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ABSTRACT 

Context. Several gamma-ray binaries show extended X-ray emission that may be associated to interactions of an outflow 
with the medium. Some of these systems are, or may be, high-mass binaries harboring young nonaccreting pulsars, in 
which the stellar and the pulsar winds collide, generating a powerful outflow that should terminate at some point in 
the ambient medium. 

Aims. This work studies the evolution and termination, as well as the related radiation, of the shocked-wind flow 
generated in high-mass binaries hosting powerful pulsars. 

Methods. A characterization, based on previous numerical work, is given for the stellar/pulsar wind interaction. Then, 
an analytical study of the further evolution of the shocked flow and its dynamical impact on the surrounding medium 
is carried out. Finally, the expected nonthermal emission from the flow termination shock, likely the dominant emitting 
region, is calculated. 

Results. The shocked wind structure, initially strongly asymmetric, becomes a quasi-spherical, supersonically expanding 
bubble, with its energy coming from the pulsar and mass from the stellar wind. This bubble eventually interacts with 
the environment on pc scales, producing a reverse and, sometimes, a forward shock. Nonthermal leptonic radiation 
can be efficient in the reverse shock. Radio emission is expected to be faint, whereas X-rays can easily reach detectable 
fluxes. Under very low magnetic fields and large nonthermal luminosities, gamma rays may also be significant. 
Conclusions. The complexity of the stellar/pulsar wind interaction is likely to be smoothed out outside the binary 
system, where the wind-mixed flow accelerates and eventually terminates in a strong reverse shock. This shock may be 
behind the extended X-rays observed in some binary systems. For very powerful pulsars, part of the unshocked pulsar 
wind may directly interact with the large-scale environment. 
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1. Introduction 

Gamma-ray binaries are compact sources that present non- 
thermal emission in the GeV and/or TeV bands and, typ- 
ically, they are also moderate-to-strong emitters in radio 
and X-rays. The radiation in gamma-ray binaries is related 
to transient or persistent outflows, which can interact either 
with themselves in the form of internal shocks or with the 
environment on different scales. The number of members of 
the class of gamma-ray binaries is growing, being at present 
around ten, including unconfirmed candidates (e.g. Albert 
et al. "MW, Aharonian ct al. 2005a, 2005b; Abdo et 

al. [2QQ£a, 200 9b^ 2009c. ,2010a, .20.10b, ,2.011; Tavani et al. 
I:j009a[ [20 09b| S abatini et al. EUTUl Williams et a l. I^OTUl 
Tam etaipnini HiU et aL l^OTTl Falcone et al. l^OTTl Corbet 
et al. lMTTll . 

Several gamma-ray binaries present hints or clear evi- 
dence of interactions between an outflow and their medium 
on large scales: Cygnus X-3, LS I -f61 303, LS 5039, 
PSR B1259— 63, and the gamma-ray emitter candidate 
Cygnus X-1 ( Heind l et al. Sanch ez-Sutil et al. 2008; 

Paredes et al.lMITl Pavlov et al l2011al Durant et al. 2011 
Martf et al. [TM51 Gallo et al. 1^0051 RusseU et al. .2007,) . 
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Cygnus X-1 and Cygnus X-3 are high-mass micro- 
quasars, in which accretion leads to jets that can interact 
with the ISM (e. g Vel azquez & Raga [MM Hein z [MI^ 
Heinz & Sunyaev 120021 Bosch-Ramon et al. 120051 Zavala 
et al. [20081 Bordas et al. [2009| Bosch- Ramon et al. [20TT|) . 
On the other hand, PSR B1259— 63 is a system formed 
by an 09.5 V star and a nonaccreting millisecond pulsar 
(Johnston et al. 119921 Negueruela et al. I201ip . Instead of 
coming from a jet, like in microquasars, the nonthermal 
emission of this source is thought to originate in the re- 
gion where the star and the pulsar winds collide (Tavani 
& Arons [Tl)97p , with its radio emission extending far away 
from the binary (Moldon et al. l2011ap . PSR B1259-63 also 
presents extended X-ray emission on scales of ^ 4—15" (i.e. 
a projected linear size of ~ 1.5 — 5 x 10^^ cm at the source 
distance of 2.3 kpc; Negueruela et al. I2011[) . with a flux 
~ 10~^^ erg cm~^ and photon index ~ 1.6 (Pavlov et 
al. I2011ap . The extended radio and X-ray emission would 
be produced in different regions but still in the shocked 
stellar-pulsar wind structure, which propagates away from 
the binary and should eventually terminate in the external 
medium. 

The nature of LS 5039 and LS I -1-61 303 is still unknown 
due to the lack of clear pulsar or accretion features, conclu- 
sive system dynamical information, or nonthermal emission 
evidence (e.g. Casares et al. I2005al Casares et al. I2005bl 
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Sarty et al [^OTTl Rea et al. l^CTTUll^CTTTl McSwain et al . 1^0111 
Ribo 120111 see also Bosch- Ramon & Khangulyan I2009P . 
In the past, several works have considered LS 5039 and 
LS I -1-61 303 good microquasar candidates because of their 
radio morphology and emission properties (e.g. Taylor et al. 
\WM Pare des et al. [10001 Massi et al. I^UOTl HOOl P aredes 
et al. '^ggi Ro mero e t al. Paredes et al. [20051 Bosch - 

Ramon et al. 120061 Massi & Kaufman Bernado I2009p . 
However, it has also been proposed that the phenomenol- 
ogy of these sources and, in particular, their milliarcsec- 
ond (mas) scale radio morphology may be more compatible 
with a nonaccreting pulsar scenario (e.g. Maraschi & Treves 
[TMli Ma rtocchia et al. [MI51 Dubus [MI51 Ch ernyakova et 
al. 120061 Sierpowska-Bartosik & Torres 120081 Sierpowska- 
Bartosik & Torres 120091 see however Romero et al. 12007P . 
Both binaries are high-mass systems, with LS 5039 harbor- 
ing an 06.5V star (Clark et al.^UUQ, a nd LS I +61 303, 
an early Be star (Hutchings & Crampton 1 1 98 1 p . 

Extended X-rays have been found in LS 5039, with an- 
gular size 9 ~ 1' — 2' (^ 2— 4x 10^® cm at 2.5 kpc; Casares et 
al. l2005aP . X-ray flux sa 9 x 10"^"* erg s~^ cm~^ and photon 
index ~ 1.9 — 3.1 (Durant et al. 1201ip , and possibly also in 
LS I -^61 303, w ith 6 - 10" (- 3 x 10^^ cm at 2 kpc; Frail k 
Hiellming [TMT|) . and X-ray flux « 2 x 10"" erg s'^ cm^^ 
(Paredes et al. [lOO'iEl) • It is noteworthy that the gamma-ray 
binari es HE SS J0632-H05 7 (Aha ronian et al. [MI7I Hinton 
et al. [20091 Skilton et al. [2U09l Falcone et al. [20T0l [20Tn 
Moldon et al. I2011bp and IF GL J1018.6-5856 (Cor bet et 
al. \Wn\ Pavlov et al. I2011bl de Ona Wilhelmi et al. [20T01) 
may host a nonaccreting pulsar as well, although a micro- 
quasar scenario cannot be discarded. 

Summarizing, among the several gamma-ray binaries 
presenting hints or evidence of interactions with the 
medium, three of them may (one for sure) host a nonac- 
creting pulsar. In addition, extended X-ray emission may 
be eventually discovered as well in HESS J0632-I-057 and 
IFGL J1018.6— 5856, another two pulsar binary candidates. 
At present, the dynamics and radiation of an outflow from a 
pulsar gamma-ray binary has not been studied far from the 
binary system. Given the recent findings of extended emis- 
sion, it seems necessary to analyze the flow evolution be- 
yond the stellar/pulsar wind region and its interaction with 
the ISM. This is the goal of this work, which is distributed 
as follows. In Sect. [H the main properties of the binary 
generated flow are described. In Sect. [31 the general prop- 
erties of the interaction between the flow and the external 
medium, the progenitor supernova remnant (SNR; young 
sources), or the interstellar medium (ISM; older sources) are 
characterized analytically, and in Sect.[4l the expected non- 
thermal emission from radio to gamma rays is computed. 
Finally, in Sect. [5l the results are discussed in the context 
of PSR B1259-63, LS I +61 303 and LS 5039 (assuming 
that the last two host a powerful pulsar), and some final 
remarks are made in Sect. [Bl All through the paper, cgs 
units are used. 



^ Rea et al. I2U1UI did not find evidence of extended X-rays 
in LS I +61 303 in longer observations than those studied in 
Paredes et al. 120071 However, in their observations the X-ray 
counts were integrated along one spatial dimension in the CCDs, 
which only allows finding an extension signal in specific direc- 
tions. 



2. Physical scenario 

2.1. The stellar-pulsar wind shock 

In Figure [1] the pulsar high-mass binary scenario is 
sketched. For a general description of the scenario on the bi- 
nary scales, see for instance Tavani & Arons (1997J. As seen 
in the figure, the stellar and the pulsar winds collide to form 
two bow-shaped standing shocks. An important parameter 
that characterizes the wind interaction region is the pulsar 
to stellar wind momentum flux ratio (e.g. Bogovalov et al. 
■2008)) : 
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where Mw-6.5 = (Mw/3 x 10^^ Mf^yr^^) is the stellar 
mass-loss rate, Wws.s = (ww/3 x 10^ cm s^ ) is the stellar 
wind velocity, and Lsd37 — (^sd/lO'^'' erg s^^) is the pulsar 
spin-down luminosity, taken here as equal to the kinetic lu- 
minosity of the pulsar wind. The stellar wind is assumed 
isotropic. The pulsar orbital velocity, of a few hundred 
km s~^, is much lower than Vvj, typically ~ 2 x 10^ cm s~^ 
for massive stars (e.g. Puis et al. 12009p . so it has been ne- 
glected in this section. However, the orbital velocity may 
play a relevant role in determining rj in some cases, e.g. 
in the presence of a stellar decretion disk, as discussed in 
Sect. 12.2.31 It is also noteworthy that the pulsar wind is 
simplified here as isotropic, although a more refined treat- 
ment should account for anisotropy (see, e.g., Bogovalov & 
Khangouhan 2002, Bogovalov et al.^UTQ. 

For reasonable 77-values < 1 (a scenario with 77 > 1 is 
briefiy discussed in Sect. 12.2?^ . the shocked wind structure 
is dominated by the stellar wind ram pressure, and points 
away from the star on the binary scales. The size of the 
wind colliding region can be characterized by the distance 
between the two-shocks contact discontinuity (CD) and the 
pulsar, 

-ftp - 7— , (2) 

which is « 7 X 10^^ cm for rj = 0.1 and a distance between 
the star and the pulsar of i?orb = 3 x 10^^ cm. The distance 
between the CD and the star is then w 2.3 x 10^^ cm. The 
CD, which separates the shocked stellar and pulsar winds, 
has an asymptotic half-opening angle (see Bogovalov et al. 
120081 and references therein) : 

a « 28.6° (4 - r;2/^) 7yi/3 . (3) 

Whereas the main contribution of momentum an mass gen- 
erally comes from the stellar wind, the energy is mostly 
provided by the pulsar wind, as seen from the pulsar to 
stellar wind luminosity ratio: 
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We notice that the wind of the pulsar prevents the latter 
to accrete from the stellar wind for a wide range of 77- values. 
As a conservative estimate, we can derive the 77- value for 
which i?p is equal to the gravitational capture radius (Bondi 
& Hoyle[I0111): 
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.(-Rorb^w)^ 

where Al is the pulsar mass, and -Roibi2.5 — (-Rorb/3 x 
10^2 cm). 
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Fig. 1. a) Sketch of the pulsar high-mass binary scenario on scales larger than the binary system. The star (S) and 
pulsar (P) shocked winds trace a spiral-shape due to the Coriolis force and the pulsar orbital motion. Eventually, both 
shocked winds mix due to Kelvin-Helmholtz instabilities, and the spiral arms merge, b) Cartoon of the bubble (case 
?7 < 1) driven by the shocked flows formed within the binary system (BS) and made of pulsar- and stellar-wind material. 
The bubble expands and accelerates, eventually terminating in a shock where its ram pressure is equal to the pressure 
of the SNR or the shocked ISM for old enough sources. This occurs at the contact discontinuity (CD). 



2.2. Shocked flow evolution 



2.2.1. Stellar wind momentum flux dominance 



After being shocked, the pulsar wind facing the star moves 
at a speed c/3, but suffers a strong reacceleration along the 
CD produced by a strong pressure gradient, related to the 
anisotropy and inhomogeneity of the stellar wind at the pul- 
sar location (Bogovalov et al. l2008|) . Because of the very low 
density in the downstream region, the pulsar wind shock is 
adiabatic (Bogovalov et al. 120081 Khangulyan et al. 2008). 
The stellar wind shock could be either adiabatic or radia- 
tive, but in either case the stellar wind velocity will be <C c. 
Therefore, the two shocked winds have a relative velocity 
~ c/3 — c, depending on how much the shocked pulsar wind 
has reaccelerated. Such a large velocity difference is likely to 
lead to the development of Kelvin-Helmholtz (KH) instabil- 
ities in the CD. This will produce entrainement of shocked 
stellar wind into the lighter and faster shocked pulsar wind. 
If the stellar wind shock is radiative, then the matter down- 
stream collapses into a dense and cool thin layer, diminish- 
ing even farther the stability of the whole shocked structure 
(e.g. Stevens et al. HM^ Pittard 12005)) . Perturbations gen- 
erated in the interface will initially move with the stellar 
shocked wind, growing on a timescale ^kh ~ l/csw, where 
I is the size of the perturbation, and the sound speed 
of the shocked stellar wind. Taking Csw ^ v^, ^ 10* cm s~^ 
and I ^ Rp = 10^^ cm, one obtains ^kh ^ 10"* s, approx- 
imately the lapse needed by the perturbations to enter in 
the nonlinear regime. This implies that instabilities, tur- 
bulence, and mixing will develop on spatial scales that are 
similar though slightly larger than those of the binary sys- 
tem. 

We now discuss the impact of the orbital motion on the 
shocked winds. To do this, we focus on three cases: the pos- 
sibly more common 77 < 1, (e.g. AI^ > 2.5 x 10~* Mq yr~^ 
for Vv, ^ 2 X 10* cm s~^ and Lgd ^ 10'^'' erg s~^); more 
briefly, the extreme case ry > 1; and under the presence of 
a stellar decretion disk. 



We analyze the dynamics of the gas flow in the coordi- 
nate system that corotates with the binary system at an 
angular velocity ft — 2tt/T, where T is the orbital period. 
The X-axis coincides with the line connecting the centers 
of the stars, the Y-axis lies in the orbital plane, and the 
Z-axis is normal to the orbital plane. In such a coordi- 
nate system, the Coriolis acceleration is Oc = —2 $7 x v. 
The stellar wind moves along the X-axis with speed ^ Ww, 
acquiring a velocity v± = 2 f2 v^, i in the Y-direction. In 
the frame of the shocked- wind flow, t — and x is 

the distance, along the X-axis, between the pulsar and a 
point farther away from the star. The dynamical pressure 
of the stellar wind in the Y-direction can be estimated as 
« p^vl « {2nt)^pW^ « (2m)2p^^ (P^ = p^vl). From 
this, the location where the pulsar wind total pressure is 
balanced by the stellar wind can be derived with 
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where pw = Pwo/(l + a;/i?orb)^, with p^o = Mw/47ri?^j.^Uw 
the stellar wind density at the pulsar location. For x ^ 
i?orb, Eq. ^ can be simplified to 



(47r)2cM„ 



and for x ^ R, 
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For typical parameters, the solution is closer to the one 
given by Eq. ([SJ, therefore the bending distance can be 
written as 



xo^Tx 10^2 Llil, v'Jl, Te M^'_',\^ cm , (9) 



r-1/2 
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where Tq ~ (T/IO^ s). This distance is much less than in the 
purely ballistic case, for which xq ^ cT — 3 x 10^^ Tg cm. 

It is noteworthy that the pulsar wind is shocked even in 
the direction opposite to the star, within a distance xq- 
Therefore, a significant fraction of the pulsar wind luminos- 
ity is reprocessed outside the binary system. Interestingly, 
this region may dominate the very high-energy output in 
pulsar gamma-ray binaries with high photon-photon ab- 
sorption (see the discussion in Bosch- Ramon et al. I2008p . 

If KH instabilities by themselves have not already mixed 
the shocked winds, isotropizing their particle, momentum, 
and energy fluxes, together with the Coriolis force they very 
likely will. Eventually, the shocked-wind flow will probably 
end up as a trans- or subsonic turbulent bubble, loaded with 
stellar wind mass and pulsar wind energy. The bubble is not 
confined because there is still a strong pressure gradient 
radially outwards. Therefore, and assuming that the whole 
process is adiabatic, any thermal, turbulent, or magnetic 
energy will eventually end up in the form of radial bulk 
motion. The role of the magnetic field should not change 
this picture much, unless the initial stellar and pulsar winds 
had a dynamically dominant magnetic field. (For the role 
of the magnetic field on binary scales; see Bogovalov et al. 
I201in The maximum velocity expansion of the bubble will 
be roughly its initial sound speed, which for a maximum 
entropy initial state (i.e. right after isotropization) is 



Dominant pulsar wind case 



Dominant stellar wind case 
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Typically, the mas radio emission in pulsar massive binaries 
will come from scales larger than xq , and thus this radiation 
should be strongly affected by the Coriolis force, as well as 
by shocked-wind mixing. The situation described in this 
section is pictured in Fig. [T] (left). 

2.2.2. Pulsar wind momentum flux dominance 

For a dominant pulsar wind, the shocked structure closes 
towards the star. In this case, the Coriolis force exerted by 
the light pulsar wind is too low and the shocked struc- 
ture moves ballistically, with a bending typical distance 
~ Uw T. Along the orbit, however, the natural bending of 
the shocked structure directly exposes it to the pulsar wind 
ram pressure, which pushes and opens the spiral to some 
degree. As for 77 < 1, for 77 > 1 the instabilities in the CD 
will also take place, and the formation of a mixing layer 
is expected to occupy part of the region of the free pul- 
sar wind. However, in the direction perpendicular to the 
orbital plane, the pulsar wind will likely propagate freely. 
The energetics of this free wind may be significant, even- 
tually terminating as a jet-like structure interacting with 
the surrounding medium (see, e.g., Bordas et al. |2009p . A 
sketch of the two cases discussed, with 77 > 1 and < 1, is 
shown in Fig. [2] 



2.2.3. The equatorial flow case 

The situation discussed so far is that of two spherical winds 
interacting with each other. However, some of the systems 
considered in this work host Be stars, which present decre- 
tion disks around the equator. These disks are thought to 
be much denser than the radiatively driven stellar wind. 




Fig. 2. Sketch of the two possibilities for the outcome of 
the interaction between the pulsar and the star winds, in 
the context of a binary system. To the right, the case for a 
very powerful pulsar, with > 1, is presented. To the left, 
the case with 77 < 1. 



The disk flow is quasi-Keplerian, and moves subsonically in 
the radial direction with velocities around ^ 1 km s~^ (e.g. 
Porter & Rivinius I2003P . Given the slow flow velocity, the 
flow speed determining the flux momentum ratio, in the 
pulsar reference frame, is near the orbital velocity of the 
pulsar. Outside the equatorial disk, the radiatively driven 
wind in Be stars is most likely polar. Such a configuration 
of the two stellar flows probably makes the geometry of the 
colliding wind region quite complex, enhancing the devel- 
opment of instabilities in the shocked-flow contact discon- 
tinuity (for simulations of this scenario on binary system 
scales, see Romero et al. 12007! and Okazaki et al. I201ip . 
The shocked disk material has a much larger density than 
the light polar wind. This implies that development of in- 
stabilities may be slower, although the sound speed of the 
shocked disk will be similar to the pulsar's orbital veloc- 
ity, the dynamical speed in this case. The high density of 
the shocked equatorial flow implies that it will likely be 
radiative, thereby potentiating instabilities and fragmenta- 
tion. Finally, the mass-loss rate of the disk is similar or 
even smaller than that of the polar wind, and thus once 
accelerated by the pulsar wind, the flow evolution on the 
largest scales should not be very different from the isotropic 
wind case. The dense and fragmented shocked disk ma- 
terial will likely introduce inhomogeneities into the larger 
scale shocked flow, although eventually the isotropization 
of the particle, momentum, and energy fluxes of the larger 
scale structure seems a reasonable guess. However, hydro- 
dynamical simulations on larger scales than the binary are 
mandatory for validating the assumption. 

2.3. Bubble large-scale evolution 



The pulsar started its life after a supernova explosion. 
For an explosion energy Esnr ^ 10^^ erg and ISM den- 
sity Pism = 1-7 X 10-23 g cm-3 {Uism = 10 cm '^), 

the SNR becomes radiative (Blinnikov et al. I1982p af- 
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ter ~ 2 X 10^ yr when achieving the radius Rc ~ 
11 (£^SNR5i/"ismi)^/^ pc, where £;snr51 = (£^snr/10^^ erg) 
and riismi = ('^-ism/lO cni~^). In this context, the binary 
system will leave the SNR after a time 



2.8i?SNR^c 
Pism 



1/5 



-7/5 . 



(11) 



i.e., isp « 8 X 10'' yr for i;p 2 x lO"^ cm s"^ 

Inside the SNR, given the high sound speed of the am- 
bient medium, the binary driven bubble cannot generate a 
forward shock. Once in the ISM, the bubble forms a slow 
forward shock. In both situations, the supersonic expand- 
ing bubble suffers a strong reverse shock to balance the 
external pressure, either from the hot SNR ejecta or the 
shocked ISM. This reverse shock is expected to have a speed 
~ Uexp in the bubble reference frame. If the pulsar is power- 
ful enough, the reverse shock can power a moderately bright 
nonthermal emitter, with a potential nonthermal luminos- 
ity Lnt as high as Lsd- This estimate is done by neglecting 
the radiation losses on the binary system scales, but as long 
as adiabatic losses are likely to dominate in the colliding 
wind region (e.g. Khangulyan et al. I2008|) . this approxi- 
mation should suffice at this stage. Thermal emission from 
the reverse shock can be discarded, given the low densities. 
Some extended thermal emission will be generated in the 
shocked ISM, but with its peak in the ultraviolet it may be 
hard to detect. Figure [T] (right) sketches the bubble termi- 
nation region. 

3. Dynamics of the bubble interacting with the 
medium 

The interaction between the bubble and its environment 
takes place in three steps. At early times, the bubble is 
still surrounded by the SNR in its adiabatic or Sedov phase 
(Sedov 1959). Later on, the SNR forward shock in the ISM 
becomes radiative, entering into the so-called snow-plow 
phase (Cox HSU Blinnikov et al. HM^ . Finally, after the 
binary leaves the SNR or the latter dissipates, the bubble 
interacts directly with the ISM. 

In the adiabatic regime, the SNR properties can be 
characterized through the SNR/ISM forward shock radius 
(Sedov [MID 
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-RSNR ~ 1.15 (£^SNR/Pi 

where tp is the pulsar age, velocity 

WSNR ~ (2/5) i?SNRAp , 

and the inner pressure is 

-£^SNR 



1/5 .2/5 



SNR 



(12) 



(13) 



(14) 



assumed here to be constant. Equilibrium between bubble 
and SNR pressures determines the location of the boundary 
between these two regions: -Rbcd ~ (3 L^d ip/107rpis„i 
and the equilibrium between bubble preshock and post- 
shock total pressures determines the location of the bubble 
reverse shock: 



^brs 



2vrFsNR,i'oxp 



(15) 



Since the hot ejecta region has a high sound velocity, the 
bubble forward shock inside the SNR quickly becomes a 
sound wave, and the bubble shape becomes spherical. 

For t > tc, the SNR enters into the radiative phase, 
so the snow-plow solution has to be used (Blinnikov et al. 
IHHS)- The SNR/ISM forward shock has now a radius of 
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a velocity of 

^^SNR = (2/7) i?SNR/ip : 

and an inner pressure of 
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where e = 0.24. As before, the bubble reverse shock radius 
can be estimated using Eq. ([T5|). 

The bubble/ISM direct interaction has a strong resem- 
blance to the interaction between a supersonic stellar wind 
and the ISM. Therefore, for this case the solution for a su- 
personic wind with continuous injection has been adopted 
(Castor et al. I1975p . This renders the ISM forward shock 
location at 



i?b « 0.76 
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,3/5 
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moving with velocity 



« (3/5) i?b/tp 



The inner pressure is now 



^b ^ Pism ^^b I 



(19) 



(20) 



(21) 



and again the preshock and postshock total pressure bal- 
ance gives the location of the bubble reverse shock: i?brs ~ 
•\/Lsd/27rPbWcxp. Since the bubble external medium has a 
low sound speed, the precise shape of the bubble is less 
clear now, although it has been assumed to be spherical, as 
discussed in Sect. 12.21 

For old enough sources, say > 10^ yr, the proper 
motion of the system can affect the ISM shock, leading to 
the formation of a bow-shaped ISM shock. In this case, the 
reverse shock is located approximately at the point where 
the bubble and ISM ram pressures become equal: 
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(22) 



Using the model just described, the relevant distances, 
velocities, and pressures can be computed for different sys- 
tem ages. Parameter values similar to those of LS 5039 
and LS I -1-61 303 have been adopted, and their values 
are presented in Table 1. The dynamical, as well as the 
radiative (see next section), results may also apply to 
PSR B1259— 63, although the comparison is less straight- 
forward, as discussed in Sect. [S] It is noteworthy that, 
as long as xq is much smaller than the largest scales of 
interaction with the medium, the size difference between 
PSR B1259— 63 and the other two sources should not in- 
troduce strong differences in the bubble propagation and 
termination. The results are shown in Figs. |3] and |H The 
breaks in reverse shock and CD distances, and pressures. 
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Parameter [unitsj 


value 


Star mass- loss rate [iVf0 yr""'^] 


3 X 10"'' 


Stellar wind velocity [cm s~^] 


2 X 10* 


Pulsar spin-down luminosity [erg s~^] 


lO^'^ 


DiMX eueigy [eigj 


-LU 


iOlVi liUIliUcl LlcllaiLy ICIli 1 


1 

±u 


NIagnetic-to-total energy density ratio 


0.01, 0.1 


Nonthermal-to-total luminosity ratio 
Star luminosity [erg s~^] 


0.01 


^039 


Stellar photon energy [eV] 


10 


Reverse shock velocity [cm s~^] 


10" 


Source distance [kpc] 


3 



Table 1. Adopted parameter values 



apparent in the figures for the SNR case at 2 x 10'' yr, 
are caused by the significant loss of internal energy (a frac- 
tion 1-e) that takes place in the adiabatic-to-radiative phase 
transition (Blinnikov et al. I1982p . 

The approach carried out here should be reasonable at 
a semi-quantitative level, but there are some caveats. First, 
these different stages of the bubble evolution are idealiza- 
tions, and mixed situations could take place. Also, the ini- 
tial conditions of the bubble were oversimplified, assuming 
that most of the complexity of the interaction on the binary 
scales is smoothed out. In particular, the importance of the 
geometry of the interacting outflows (e.g. polar winds, disk, 
etc.) is to be studied in detail. On larger scales, the shocked 
bubble, SNR and ISM media were approximated as homo- 
geneous static gas. Rayleigh- Taylor instabilities in the CD 
between the denser shocked ISM shell, and its inner re- 
gion, were neglected, along with anisotropics in the ISM. 
To account for all this properly, magnetohydrodynamical 
simulations of the bubble formation, evolution and termi- 
nation are planned. We also neglected the pulsar spin-down 
luminosity evolution with time, which may imply an under- 
estimate of the total bubble injected energy by a factor of 
a few. 



4. Nonthermal radiation estimates 

The reverse shock is a good candidate for nonthermal emis- 
sion, since it is fast and strong, and basically all the bubble 
energy goes through it. The forward shock in the ambi- 
ent medium, on the other hand, is either much slower and 
also less energetic (bubble/ISM case), or just a sound wave 
(SNR case). 

4.1. Emitter characterization 

To compute the particle acceleration rate in the reverse 
shock, we adopted the expression for a non-relativistic hy- 
drodynamical shock in the test particle and Bohm diffu- 
sion approximations: E = (l/27r)(wexp/c)^ qBc (e.g. Drury 
I1983p . where B is the magnetic field and q the particle 
charge. The same B and diffusion coefficient values were 
assumed in both sides of the shock. The resulting energy 
distribution of the particles injected in the emitter depends 
on energy as oc E^'^. The _B-energy density was taken equal 
to ^SNR|b7 with = 0.1 and 0.01. Smaller values seem 
less plausible, because some magnetic field is expected to 



a) 




Ie+I9 - 

( 

- ; 
I - -'' 

le+18 - 



SNR forward shock 

- ■ Biibble/SNR CD 
. . ■ Bubble reverse shock 



le+17 I • 1 — I'll 1 • • • — I'll 

10000 le+05 

b) 

le+20 u • • 




a: 



le+18 - 



I I , , , , , III 

10000 le+05 



Fig. 3. a) Evolution with time of the radn of the SNR/ISM 
forward shock (black solid line), SNR/bubble contact dis- 
continuity (red dashed line), and bubble reverse shock (blue 
dotted line). The CD radius becomes larger than the SNR 
shock radius, which indicates that the model adopted here 
does not apply after a few x 10^ yr. b) Evolution with time 
of the radii of the bubble/ISM forward (black solid line) 
and reverse shocks (blue dotted line). 



be transported from the binary scales, but cannot be dis- 
carded. The maximum energies are the minimum value be- 
tween the synchrotron-limited case, i^max ~ \/ Q c/ag B, 
with Os = 1.6 X lO^'^, and the diffusive-escape one, -Emax ~ 
-^3/477 (fb/c) q B i?brs- Although the stellar photon energy 
density, the dominant radiation field, is higher than the B- 
energy density at i?brs, the Klein-Nishina (KN) effect makes 
IC losses less important than synchrotron (or escape) ones 
at i?max- For a wide source age range, in the case of electrons 
Emax is ~ 100 TeV, but determined by synchrotron cooling 
for ^B = 0.1, and by diffusive escape for ^b = 0.01. For 
protons, the maximum energies are 300 TeV for ^b — 0.1 
and 100 TeV for ^b = 0.01, both limited by diffusive escape. 
The differences between the maximum energies in the SNR 
and the bubble/ISM case are small; in particular, under 
diffusive escape, they are equal. Synchrotron-limited max- 
imum energies are cx but diffusive -escape ones are 
constant, because B x _Rbrs is also constant. 
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a) 
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Fig. 4. a) Evolution with time of the SNR/ (black solid 
line) and bubble/ISM forward shock velocities (red dashed 
line), b) Evolution with time of the SNR/ (black solid line) 
and bubble/ISM pressures (red dashed line), c) Evolution 
with time of the SNR/ (black solid/blue dashed line) and 
bubble/ISM magnetic fields (red dotted/green dot-dashed 
line). 

The conditions downstream of the bubble reverse shock, 
in particular the low densities, render hadronic processes in- 
efficient. Therefore, we have focused on electrons emitting 



via synchrotron and inverse Compton (IC) scattering. The 
most energetic protons may diffuse away from the reverse 
shock to reach the shocked ISM, but to be efficient, this 
mechanism requires target densities to be much higher than 
those considered here. To compute the synchrotron and IC 
emission (e.g. Blumenthal & Gould [TOTOI . we modeled the 
emitter as only one zone, taking homogeneous B- and ra- 
diation fields. This approach is valid for the synchrotron 
calculations as long as the flow is subsonic, and is the 
same everywhere, which is admittedly a rough approxima- 
tion in an ideal MHD flow. The role of diffusion in electron 
transport is negligible in the present context unless diffu- 
sion coefficients are much higher than Bohm. Owing to fast 
cooling, synchrotron X-rays are radiated close to the re- 
verse shock region, whereas the synchrotron radio emission 
comes from a more distant location due to advection. The 
particle flux conservation means that the advection speed 
downstream of the reverse shock is roughly oc 1/ B? , where 
R is the distance to the reverse shock. This makes particles 
accumulate at 

R ~ i?brs (1 + 3Wexptcool/4i?brs)'/' ^ 10 i?brs , (23) 

with tcooi the relevant cooling time. Most of the IC radia- 
tion comes also from a distance ^ i? to the star. This is not 
true for very high-energy IC photons, but for the adopted 
B-values and due to the KN effect, their fluxes are minor so 
were neglected here (see however Sect. [5]) . Given the quasi- 
spherical shape of the emitter, the target photon field for 
IC has been taken as isotropic and monoenergetic, given 
its narrow band. The stellar luminosity and photon energy 
have been fixed to 10"^^ erg s~^ and 10 eV. The luminosity 
injected in the form of nonthermal electrons at the bubble 
reverse shock was taken as 1% of the pulsar wind lumi- 
nosity: I/nt = 0.01 Lsd- Although Lnt is hard to determine, 
changes in this parameter only affect the nonthermal fluxes 
linearly. Adiabatic cooling has been included in the calcu- 
lations, with the cooling timescale taken as i?SNR|b/'i'SNR|b- 

4.2. Spectral energy distributions and lightcurves 

After characterizing the emitter as a region with roughly 
homogeneous properties, the maximum particle energies, 
and the dominant cooling processes (adiabatic and radia- 
tive), it is easy to calculate the electron evolution. For that, 
a constant particle injection is assumed, and electrons are 
evolved up to the source age of interest. Once the particle 
energy distribution is known, the synchrotron and IC spec- 
tral energy distributions, as well as the specific and inte- 
grated fluxes, are calculated by convolving the synchrotron 
and IC speciflc power per electron by the particle energy 
distribution. 

The results of the radiation calculations are shown in 
Figs.[5]and[6l Figure[5]shows the computed synchrotron and 
IC spectral energy distributions for a source age of 10* yr 
(SNR) and 3 x 10* yr (bubble/ISM), with the remaining 
parameter values given in Table 1. Figure |6] presents the 
time evolution of the surface brightness at 5 GHz, and the 1- 
10 keV and 0.1-10 GeV bolometric fluxes. The radio surface 
brightness was computed by just dividing the specific flux 
by the source solid angle in 1" x 1" units at 3 kpc, taking 
R as the source size. We point out that our aim is just to 
provide with a crude estimate of the expected radiation. 
More detailed calculations of the nonthermal emission will 
be presented elsewhere. 
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As seen in the figures, the break of the adiabatic- 
to-radiative phase transition appears in ah the SNR 
hghtcurves, although is more apparent for the radio surface 
brightness and IC fluxes, which are strongly affected by the 
jump in R. The general evolution of R produces the long- 
term strong decay of the SNR radio and GeV lightcurves, 
and a smoother decay of the bubble/ISM radio lightcurve, 
whereas the bubble/ISM GeV lightcurve remains more or 
less constant. The differences are partially caused by the 
evolution of the relative importance of the different radia- 
tion channels and the adiabatic one, which is not the same 
in the SNR and the bubble/ISM cases. The synchrotron and 
IC emissivities have also different time dependencies, the 
former depending on as well. For = 0.1, X-ray emit- 
ting electrons radiate all their energy through synchrotron, 
rendering a flat X-ray lightcurve and photon indexes ^ 2. 
For ^B = 0.01, with the electron maximum energies limited 
by diffusive-escape, the maximum energy of synchrotron 
photons {oc B E^) decreases with time, and it quenches the 
X-ray flux for tp ^ 10^ yr. This also explains the steepening 
of the synchrotron X-ray spectrum for the bubble/ISM case 
seen in Fig.[5l The SNR X-ray spectrum does not show this 
feature because it has been computed for t = 10'* yr < tc- 
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Fig. 5. Spectral energy distributions of the synchrotron 
(thick lines) and IC (thin lines) emission computed for the 
SNR (10"* yr) and the bubble/ISM (3 x 10"* yr) cases, and 
Cb =0.01, 0.1. 



l-IOkeV (SNR)5j=0.1 

l-IO keV (Bubble/ISM) 5^=0.1 

l-I0keV(SNR)5|j=0.01 

1-10 keV (Bubble/ISM) 5j=0.01 



0.1-10 GeV (SNR) 5j=0.I 

_. 0.1-10 GeV (Bubble/ISM) 5^=0.1 

. . 0.1-10 GeV (SNR) 5|j=0.01 

. - 0.1-10 GeV (Bubble/ISM) E =0.01 




5. Discussion 

The sizes of the bubble reverse shock and its downstream 
region give an idea of the source angular size at differ- 
ent wavelengths. At 3 kpc, the hard X-ray emitter would 
have e ^ 10" when stiU inside the SNR {t 10^ yr), and 
- 20" - 1' in the bubble/ISM case {t « 10* - 10^ yr). 
These 9 values were calculated for the pism, -E-snr, and Lgd 
values given in Table 1, although we need to recall their 
weak inter-dependence: 9 (x (-Esnr | -^sd/pism)^^^- Extended 
X-ray fluxes could be detectable by an instrument like 
Chandra up to tp ~ 10^ yr. It is worth noting that, if 



Fig. 6. a) Lightcurve of the radio (5 GHz) surface bright- 
ness, b) Lightcurve of X-ray flux in the range 1-10 keV. c) 
Lightcurve of the IC flux in the range 0.1-10 GeV. Both 
the SNR and the bubble/ISM cases are shown in the age 
range ip = 3 x (10^ - 10^) yr and for ^b = 0.01, 0.1. 



synchrotron X-ray energies are detected from the reverse 
shock, the diffusion coefficient should be close to Bohm. 
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Radio emission is possibly detectable for the SNR case 
with = 0.1, it may be hard to be found for the case — 
0.01, and its detection seems discarded for the bubble/ISM 
case for any ^B-value. 

The GeV fluxes seem too low to be detectable by any 
current instrument, and the chances are even lower in TeV 
due to the Klein-Nishina effect and dominant synchrotron 
cooling. However, it may still be possible to get higher 
IC GeV and TeV fluxes increasing ten times Lnt, and re- 
ducing strongly ^b not to overcome the radio/X-ray con- 
straints. This would lead to GeV- TeV fluxes that could be 
detectable by present (GeV: Fermi, AGILE; TeV: HESS, 
MAGIC, VERITAS) and/or forthcoming instruments (e.g. 
TeV: CTA). Given the spectral break at - 10 GeV intro- 
duced by adiabatic losses, the bubble reverse shock would 
be a good target for a ~ 10 GeV-threshold Cherenkov in- 
strument. 

We now briefly discuss whether the evidence or hints 
of extended X-rays from LS 5039, LS I -1-61 303, and 
PSR B1259— 63 can be understood in the scenario posed 
here. 

5.1. LS 5039 

LS 5039 is likely a young source with an age range (4 — 
10) X 10^ yr, and it has already abandoned its SNR (e.g. 
Ribo et al. 2002; Moldon et al. 1201 lc| . The angular size 
at 2.5 kpc for the bubble/ISM case is consistent with the 
observed angular size of the emitter, ~ 1' — 2', and the 
fluxes can be explained in the framework given here. The 
photon index at distances ^ 1' is ^ 1.9 ± 0.7, consistent 
with the SED shown in Fig. [SI and then softens farther out 
reaching - 3.1 ± 0.5 at ~ 2' (Durant et al. l^OTTI) . This can 
be explained by the cooling of the highest energy electrons 
close to the reverse shock, thereby depleting the highest 
energy part of photons and softening the spectrum farther 
from the shock. The adopted riism = 10 cm~'^ is similar 
to the values derived for the surroundings of the LS 5039 
SNR candidate (Ribo et al. I2002[) . and the pulsar Lsd is 
probably ~ lO'^^ erg s~^, which is required to explain the 
GeV luminosity of the source (see Sect. 4 in Zabalza et 
al l2011a|) . The morphology of the extended radio emission 
found by Durant et al. pOlip is not completely spherical, 
which may be explained by the only partial isotropization 
of the momentum flux in the inner regions of the bubble 
(see Sect. 12.2)1 . A scenario with 77 > 1 could be also behind 
the anisotropy (see Sect. 12. 2T^ . but it seems less likely. 

5.2. LS I +61 303 

Assuming that the hints of extended X-rays in LS I +61 303 
are real, we see that the angular size, ~ 10" , is compatible 
with the bubble/SNR scenario for tp < 10^ yr, but with 
a low nonthermal efficiency or a very low magnetic field, 
given the weak X-ray flux (Paredes et al. I2007P and the 
radio nondetection on those scales (Marti et al. I1998P . 

The SNR progenitor of LS I -1-61 303 has not yet been 
found. However, extended 5 GHz emission has been de- 
tected centered on the location of LS I -|-61 303, with a typ- 
ical size of 6-8 arcminutes and fluxes of tens of mJy (Marti 
et al. 119981 Mufioz-Arjonilla et al. .2009 ) . It has been argued 
that the lack of extended X-rays on the same scales, and 
the low surface radio brightness (in a nonthermal scenario) , 



may go against an SNR interpretation (Marti et al. 119981 
Mufioz-Arjonilla et al. 12009^ . However, the radio emission 
surrounding LS I -1-61 303 may be mostly thermal (Muhoz- 
Arjonilla et al. l2009p . coming from an SNR/ISM shock close 
to its radiative phase. If the nism values in the region of 
LS I +61 303 were a bit higher than those adopted here, an 
SNR remnant with energy i?sNR ~ 10^^ erg would enter its 
radiative phase after 10* yr, with a radius 10^^ cm, a 
few arcminutes at 2 kpc. Thermal X-rays would not be ex- 
pected in that case, since the shocked ISM material would 
be too cold. The slow proper motion of the source (Dhawan 
et al. I2006P is compatible with LS I +61 303 being at the 
center of the SNR (see however Marti et al. JJi98 for image 
deformation). As in LS 5039, the Lsd-value in LS I +61 303 
is expected to be ~ 10'^'' erg s~^ because of the hi gh GeV 
luminosity of the source (see Sect. 5 in Zabalza et al l2011bp . 

5.3. PSR B1259-63 

PSR B1259-63 has an age of tp « 3 x 10^ yr (Wex et al. 
ri99^ . and is likely now interacting directly with the ISM. 
With this age, and the pulsar spin-down luminosity Lgd ~ 
8 X 10^5 erg s'^ (Manchester et al. [TM5)) . Eq. (HH) yields 
i?brs ~ 2 X 10^^ cm (~ 1'), assuming nism = 10 cm^^. This 
is about one arcminute at 2.3 kpc and much larger than the 
observed 6* ~ 4" -15" (Pavlov et al. l2011ap . Nevertheless, in 
PSR B1259— 63 it seems likely that the proper motion has 
already bow-shaped the ISM shock. In that case, for Tiism = 
10 cm~'^ and fp ~ 10^ cm s~^, the reverse shock would be 
located at a distance of ^ 3 x 10^^ cm from the system, 
or ~ 10" neglecting projection effects. Unfortunately, the 
large errors of the proper motion measurements (Zacharias 
et al. l2009p and the low statistics of the X-ray data make it 
difficult to compare the X-ray extension and proper motion 
directions. 

There is another possibility behind the extended X-rays 
found in PSR B1259— 63. The most significant detection 
comes from ~ 4" , or ~ 10^^ cm in linear (projected) size, 
which is ^ VcxpT- Although the flow bending distance, xq, 
should be less than that, the X-ray emission may be related 
to the asymmetric shock produced by the Coriolis force, 
and/or to spiral-arm merging. This radiation could not be 
resolved in the case of LS I +61 303 and LS 5039, which 
are much more compact than PSR B1259— 63. 



6. Final remarks 

Gamma-ray binaries hosting powerful pulsars can produce 
supersonic flows that originate in the interaction of a pulsar 
wind with the wind of the companion, which may present 
strong anisotropics and inhomogeneities. The likely mixing 
of these flows on larger scales than the binary will render 
an expanding, rather isotropic, supersonic bubble, which 
will eventually terminate in the surrounding medium. The 
interaction of the supersonic bubble with the surrounding 
medium can take place in different ways depending on the 
age of the pulsar. For young sources, the environment will 
be the hot SNR ejecta, whereas it will be the shocked ISM 
for older objects. Although this interaction is expected to 
be rather symmetric, for old enough sources, possibly like 
PSR B1259— 63, the proper motion can bow-shape the in- 
teraction structure with the ISM. 
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The radiation from the shocked bubble interacting 
with the environment may explain the observed extended 
X-ray emission from LS 5039, and possibly also from 
LS I +61 303 (if real). For PSR B1259-63, the found ex- 
tended X-rays could come from inner regions of the bubble, 
triggered perhaps by Coriolis force shocks or some other 
type of dissipation mechanism. Extended X-ray emission 
may also eventually be detected in HESS J0632-t-057 and 
IFGL J1018.6— 5856, which could also originate in shocked 
wind outflows. The shape of the interaction region might 
allow the distinction of the driving flow, i.e. a jet or a pulsar 
wind, but as shown in Sect. I^T^TTl it may be not straightfor- 
ward. 

The complexity of the flow evolution on the different 
relevant scales makes a proper analytical characterization 
difficult, and thus makes magnetohydrodynamical simula- 
tions important. The geometry, level of anisotropy, veloc- 
ity and density of the stellar flow requires careful study, 
in particular for binaries hosting stars with an equatorial 
disk, like PSR B1259-63 and LS I +61 303. Anisotropy in 
the pulsar wind has been also neglected here, but may play 
some role, at least up to intermediate scales. 
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